#LIBRERIAS
library("skimr")
library(fdth)
library(naniar)
library(BHH2)
library(agricolae)

Evaluación del recurso eólico

Proyecto de Astroestadística, Chávez Bartolo Julieta Janeth

INTRODUCCIÓN

& Problema comercial

Actualmente resolver ecuaciones de fluidos para lugares abiertos parece ser una tarea imposible. Por ello la ciencia de datos nos brinda una herramienta poderosisíma que está dando grandes resultados en la vida moderna, si bien, no resolveré un método en especifíco con este trabajo mi objetivo será mostrar lo que es posible lograr con un conjunto de datos de un fluido como es el aire.

Hoy día la producción de energía eólica ya parece ser un negocio viable, pues sus costos han disminuido considerablemente gracias al desarrollo tecnológico y a la implementación de proyectos para la evaluación del recurso. Como vemos en la Figura 1.

texto_alternativo

Sin embargo el costo de un proyecto de instalación de turbina eólica sigue siendo alto en cualquier sentido, desde la máquina hasta la mano de obra. Por ello es necesaria una correcta evaluación del recurso eólico en la zona antes de aventurarnos a instalar una turbina ya que de no hacerlo las pérdidas económicas podrían ser devastadoras e irreversibles.En la figura 2 podemos ver los porcentajes de costos como de inversión inicial.

texto_alternativo

Curva de Potencias

Esta es la relación, proporcionada por el fabricante, que existe entre la velocidad de viento y la potencia que se podría producir. Claro que esto en la vida real no es tan simple, pues el tipo de flujo de viento dependerá enormemente. Así que es una mera simplificación. El estudio de tipo de flujo es un estudio a parte. En este trabajo estaré considerando que el tipo de flujo es laminar durante el tiempo de medición.

El aerogenerador que se elija deberá tener una curva de potencias que encaje con la distribución de velocidades de viento del citio en el que se desee poner el aerogenerador. Ya que por ejemplo, si se omite este paso es posible que gastemos al rededor \(20\) millones de pesos en la turbina, pero sea inservible para el lugar porque podría ser que las velocidades del sitio nunca sean suficientes para la velocidad de arranque o bien es posible que sean demasiado fuertes como para romperlo, forzando al aerogenerador a estar detenido la mayor parte del tiempo como mecanismo de seguridad. Este es el ejemplo perfecto de cuando se tiene un gran recurso eólico, pero no se puede aprobechar a causa de un mal estudio. En la siguiente figura podrémos ver una ejemplificación general de la curva de potencias de un aerogenerador.

texto_alternativo

Datos

Los datos me fueron proporcionados por el Dr.Oswaldo Rodriguez del instituto de energías renovables,Morelos, UNAM. Estos datos provienen de un aerogenerador ubicado en una granja eólica en Oaxaca.

Las mediciones se realizaron durante el invierno, en el mes de noviembre de 2017. Y se registraron cada 10 minutos, eso nos da 30X24X6= 4320 mediciones. Al registrarse datos en un intervalo de 10 segundos, se registran cantidades de las cuales se guarda la máxima, la mínica, el promedio y la desviación estándar. Donde las columnas son:

YYYY: Representa el año, 2017. El tipo de variable es cuantitativa y su escala es de intervalo.

DD: El día. El tipo de variable es cualitativa y su escala nominal. Se usan códigos numéricos.

**MM:* El mes, noviembre denotado por 11. El tipo de variable es cualitativa y su escala es ordinal.

hh: Denota la hora entera. Va de [00,23] en unidades enteras. El tipo de variable es cuantitativa y su escala es de intervalo.

mm: Denota los minutos. Va de [00,50] en unidades de 10 en 10.El tipo de variable es cuantitativa y su escala es de intervalo. Variable discreta.

WS_Xm_mean: Velocidad del viento a la altura X en metros. Hay registros para X=80m, X=60m,X=40m y X=20m.A la altura de 80m se tienen 2 datos, de la torre A y la torre B. El tipo de variable es cuantitativa y su escala es de razon. Variable continua.

WD_X-2m_mean: La potencia producida en kW/hora, notamos que es X-2 ya que no necesariamente el medidor de velocidad está al centro de la turbina eólica. En el modelo Vestas de aerogenerador está 2 metros abajo. Solo hay registros para X=60m y X=80m.El tipo de variable es cuantitativa y su escala es de razon. Variable continua.

temp_Ym_mean: la temperatura en grados centígrados Para estas cantidades notemos que están medidas a Y=15m, Y=40m y Y=80m. De modo que solo había un único medidor a esta altura. El tipo de variable es cuantitativa y su escala es de razon. Variable continua.

RH_15m_mean: la presión en hectopascales. (1 hPa= 100 Pa). La humedad relativa El tipo de variable es cuantitativa y su escala es de razon. Variable continua.

P_15m_mean: Representa la irradiancia. La radiación solar medida en cada una de las estaciones meteorológicas es ofrecida en unidades de potencia y está en vatios por metro cuadrado (W/m²) El tipo de variable es cuantitativa y su escala es de razon. Variable continua.

Rad_80m_mean: intenté explorar esta columna para darme una idea de lo que pudiera ser, pero no tuve éxito ya que son puros ceros. Al investigar esta medida Rad, en energía eólica representa los giros por segundo de las aspas. De modo que es un dato que no usaría pues evidentemente tiene problemas de registro. De modo que será eliminada.

DENSIDAD_15m_mean: la densidad del viento [Kg/m^3]. El tipo de variable es cuantitativa y su escala es de razon. Variable continua.

BATERIA_mean: realmente desconozco la función de esta columna. Será eliminada.

1 Ańalisis Exploratorio de Datos

Cargemos el conjunto de datos

data <- read.table(file = "M06_201711.txt", header = TRUE, fill=TRUE)
head(data)
##     YYYY  MM  DD  hh  mm WS_80mA_mean WS_80mA_min WS_80mA_max WS_80mA_stdev
## 1 Height  80  80  80  80           80          80          80            80
## 2  Units m/s m/s m/s m/s          m/s         m/s         m/s           m/s
## 3   2017  11  01  00  00   7.9980E+00  5.7843E+00  9.8262E+00    8.1147E-01
## 4   2017  11  01  00  10   7.5565E+00  5.4734E+00  8.8934E+00    7.1261E-01
## 5   2017  11  01  00  20   7.5254E+00  5.7843E+00  8.8934E+00    6.3488E-01
## 6   2017  11  01  00  30   6.9534E+00  5.1625E+00  7.9607E+00    6.1436E-01
##   WS_80mB_mean WS_80mB_min WS_80mB_max WS_80mB_stdev WS_60m_mean WS_60m_min
## 1           60          60          60            60          40         40
## 2          m/s         m/s         m/s           m/s         m/s        m/s
## 3   8.0441E+00  6.1043E+00  9.5239E+00    7.9147E-01  7.0153E+00 4.8624E+00
## 4   7.5467E+00  5.7934E+00  9.2130E+00    6.8951E-01  6.7975E+00 5.4846E+00
## 5   7.4908E+00  5.7934E+00  8.9021E+00    6.3107E-01  6.6171E+00 5.1735E+00
## 6   6.9436E+00  5.1717E+00  7.9695E+00    5.9687E-01  6.1193E+00 4.2402E+00
##   WS_60m_max WS_60m_stdev WS_40m_mean WS_40m_min WS_40m_max WS_40m_stdev
## 1         40           40          20         20         20           20
## 2        m/s          m/s         m/s        m/s        m/s          m/s
## 3 8.9068E+00   9.1963E-01  6.0012E+00 3.9260E+00 8.2752E+00   8.9967E-01
## 4 8.2846E+00   7.2674E-01  5.6346E+00 3.9260E+00 7.0326E+00   6.6916E-01
## 5 8.5957E+00   7.4479E-01  5.7340E+00 3.9260E+00 7.9646E+00   7.6857E-01
## 6 7.6624E+00   6.7261E-01  5.2929E+00 3.6153E+00 7.0326E+00   7.0582E-01
##   WS_20m_mean WS_20m_min WS_20m_max WS_20m_stdev WD_78m_mean WD_78m_min
## 1          78         78         78           78          58         58
## 2           °          °          °            °           °          °
## 3  4.9511E+00 2.6832E+00 7.0361E+00   8.0715E-01  3.1790E+02        NaN
## 4  4.7894E+00 2.9941E+00 6.7252E+00   7.7481E-01  3.2390E+02        NaN
## 5  4.7602E+00 3.3051E+00 7.3470E+00   7.1574E-01  3.2570E+02        NaN
## 6  4.5873E+00 3.3051E+00 6.1034E+00   6.1686E-01  3.2710E+02        NaN
##   WD_78m_max WD_78m_stdev WD_58m_mean WD_58m_min WD_58m_max WD_58m_stdev
## 1         58           58          80         80         80           80
## 2          °            °          °C         °C         °C           °C
## 3        NaN   4.5610E+00  3.1290E+02        NaN        NaN   6.0170E+00
## 4        NaN   4.3810E+00  3.1970E+02        NaN        NaN   5.5730E+00
## 5        NaN   3.7020E+00  3.2090E+02        NaN        NaN   4.6730E+00
## 6        NaN   3.7040E+00  3.2150E+02        NaN        NaN   4.3770E+00
##   temp_80m_mean temp_80m_min temp_80m_max temp_80m_stdev temp_40m_mean
## 1            40           40           40             40            15
## 2            °C           °C           °C             °C            °C
## 3    2.4260E+01   2.4160E+01   2.4370E+01     3.6000E-02    2.4170E+01
## 4    2.3990E+01   2.3960E+01   2.4230E+01     4.2000E-02    2.3900E+01
## 5    2.4160E+01   2.4030E+01   2.4230E+01     3.4000E-02    2.4090E+01
## 6    2.4040E+01   2.3960E+01   2.4160E+01     3.6000E-02    2.4000E+01
##   temp_40m_min temp_40m_max temp_40m_stdev temp_15m_mean temp_15m_min
## 1           15           15             15            15           15
## 2           °C           °C             °C             %            %
## 3   2.4090E+01   2.4230E+01     3.5000E-02    2.4020E+01   2.3960E+01
## 4   2.3890E+01   2.4160E+01     3.3000E-02    2.3760E+01   2.3690E+01
## 5   2.4030E+01   2.4160E+01     2.8000E-02    2.3950E+01   2.3890E+01
## 6   2.3890E+01   2.4090E+01     4.4000E-02    2.3870E+01   2.3750E+01
##   temp_15m_max temp_15m_stdev RH_15m_mean RH_15m_min RH_15m_max RH_15m_stdev
## 1           15             15          15         15         15           15
## 2            %              %         hPa        hPa        hPa          hPa
## 3   2.4090E+01     2.8000E-02  8.4100E+01 8.4000E+01 8.4300E+01   6.2000E-02
## 4   2.4030E+01     2.8000E-02  8.4200E+01 8.3900E+01 8.4500E+01   1.3000E-01
## 5   2.4030E+01     2.8000E-02  8.4700E+01 8.4400E+01 8.5000E+01   1.6600E-01
## 6   2.3960E+01     3.8000E-02  8.5200E+01 8.4900E+01 8.5500E+01   1.5300E-01
##   P_15m_mean P_15m_min P_15m_max P_15m_stdev Rad_80m_mean Rad_80m_min
## 1         80        80        80          80           15          15
## 2       W/m2      W/m2      W/m2        W/m2       kg/m^3      kg/m^3
## 3 1.0000E+03       NaN       NaN  9.1000E-02   0.0000E+00  0.0000E+00
## 4 1.0000E+03       NaN       NaN  6.9000E-02   0.0000E+00  0.0000E+00
## 5 1.0000E+03       NaN       NaN  9.1000E-02   0.0000E+00  0.0000E+00
## 6 1.0000E+03       NaN       NaN  8.1000E-02   0.0000E+00  0.0000E+00
##   Rad_80m_max Rad_80m_stdev DENSIDAD_15m_mean DENSIDAD_15m_min DENSIDAD_15m_max
## 1          15            15                15               15               15
## 2      kg/m^3        kg/m^3                 V                V                V
## 3  0.0000E+00    0.0000E+00        1.9830E+00       1.9740E+00       1.9920E+00
## 4  0.0000E+00    0.0000E+00        2.0080E+00       1.9830E+00       2.0140E+00
## 5  0.0000E+00    0.0000E+00        2.0040E+00       1.9920E+00       2.0170E+00
## 6  0.0000E+00    0.0000E+00        2.0240E+00       2.0080E+00       2.0390E+00
##   DENSIDAD_15m_stdev BATERIA_mean BATERIA_min BATERIA_max BATERIA_stdev
## 1                 15           NA          NA          NA            NA
## 2                  V           NA          NA          NA            NA
## 3         3.0000E-03        12.85         NaN         NaN           NaN
## 4         4.0000E-03        12.69         NaN         NaN           NaN
## 5         6.0000E-03        12.83         NaN         NaN           NaN
## 6         6.0000E-03        12.84         NaN         NaN           NaN

Eliminemos las columnas llenas de NANś.

Para no estar cargando con los datos, quedémonos solamente con los registros promedio de cada variable. Esto lo hago bajo la justificación de que es posible ignorar el error de cada medición porque este error, denotado por la desviación estándar, será mucho menor la desviación estándar del conjunto de datos. Támbién se eliminarán las variables con problemas de registo.

data = data[-1,]
data = data[-1,]
data$WS_20m_min = NULL
data$WS_20m_max = NULL
data$WS_20m_stdev = NULL
data$WS_40m_min = NULL
data$WS_40m_max = NULL
data$WS_40m_stdev = NULL
data$WS_60m_min = NULL
data$WS_60m_max = NULL
data$WS_60m_stdev = NULL
data$WS_80m_min = NULL
data$WS_80m_max = NULL
data$WS_80m_stdev = NULL
data$WS_80mA_min = NULL
data$WS_80mA_max = NULL
data$WS_80mA_stdev = NULL
data$WS_80mB_min = NULL
data$WS_80mB_max = NULL
data$WS_80mB_stdev = NULL

data$WD_78m_min = NULL
data$WD_78m_max = NULL
data$WD_78m_stdev = NULL
data$WD_58m_min = NULL
data$WD_58m_max = NULL
data$WD_58m_stdev = NULL
data$P_15m_min = NULL
data$P_15m_max = NULL
data$BATERIA_min = NULL
data$BATERIA_mean = NULL
data$BATERIA_max = NULL
data$BATERIA_stdev = NULL
data$DENSIDAD_15m_max = NULL
data$DENSIDAD_15m_min = NULL
data$DENSIDAD_15m_stdev = NULL
data$Rad_80m_min = NULL
data$Rad_80m_max = NULL
data$Rad_80m_stdev = NULL
data$Rad_80m_mean = NULL
data$P_15m_min = NULL
data$P_15m_max = NULL
data$P_15m_stdev = NULL
data$RH_15m_stdev = NULL
data$RH_15m_min = NULL
data$RH_15m_max = NULL
data$temp_15m_stdev = NULL
data$temp_15m_min = NULL
data$temp_15m_max = NULL
data$temp_40m_stdev = NULL
data$temp_40m_min = NULL
data$temp_80m_max = NULL
data$temp_80m_stdev = NULL
data$temp_80m_min = NULL
data$temp_80m_max = NULL
data$YYYY = NULL
head(data,3)
##   MM DD hh mm WS_80mA_mean WS_80mB_mean WS_60m_mean WS_40m_mean WS_20m_mean
## 3 11 01 00 00   7.9980E+00   8.0441E+00  7.0153E+00  6.0012E+00  4.9511E+00
## 4 11 01 00 10   7.5565E+00   7.5467E+00  6.7975E+00  5.6346E+00  4.7894E+00
## 5 11 01 00 20   7.5254E+00   7.4908E+00  6.6171E+00  5.7340E+00  4.7602E+00
##   WD_78m_mean WD_58m_mean temp_80m_mean temp_40m_mean temp_40m_max
## 3  3.1790E+02  3.1290E+02    2.4260E+01    2.4170E+01   2.4230E+01
## 4  3.2390E+02  3.1970E+02    2.3990E+01    2.3900E+01   2.4160E+01
## 5  3.2570E+02  3.2090E+02    2.4160E+01    2.4090E+01   2.4160E+01
##   temp_15m_mean RH_15m_mean P_15m_mean DENSIDAD_15m_mean
## 3    2.4020E+01  8.4100E+01 1.0000E+03        1.9830E+00
## 4    2.3760E+01  8.4200E+01 1.0000E+03        2.0080E+00
## 5    2.3950E+01  8.4700E+01 1.0000E+03        2.0040E+00

Antes de hacer operaciones cambiemos el tipo de dato.

for (i in 1:18) {
    print(class(data[,i] ))
}
## [1] "character"
## [1] "character"
## [1] "character"
## [1] "character"
## [1] "character"
## [1] "character"
## [1] "character"
## [1] "character"
## [1] "character"
## [1] "character"
## [1] "character"
## [1] "character"
## [1] "character"
## [1] "character"
## [1] "character"
## [1] "character"
## [1] "character"
## [1] "character"

Cada columna tiene tipo de dato character, de modo que lo cambiaré a numeric para facilitar su manejo.

Hagamos un ciclo for en el que convierta cada columna en numeric.

for (i in 1:18) {
    data[,i] <- as.numeric(data[,i] )
}

Veamos que efectivamente se hizo el cambio de variable.

for (i in 1:18) {
    print(class(data[,i] ))
}
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
attach(data) # con el attach ya puedo accesar a los nombres de las variables
library("skimr")
skimr::skim(data)
Data summary
Name data
Number of rows 4321
Number of columns 18
_______________________
Column type frequency:
numeric 18
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
MM 0 1 11.00 0.02 11.00 11.00 11.00 11.00 12.00 ▇▁▁▁▁
DD 0 1 15.50 8.66 1.00 8.00 15.00 23.00 30.00 ▇▇▇▇▇
hh 0 1 11.50 6.92 0.00 5.00 11.00 17.00 23.00 ▇▇▆▇▇
mm 0 1 24.99 17.08 0.00 10.00 20.00 40.00 50.00 ▇▃▃▃▃
WS_80mA_mean 0 1 6.64 3.12 0.19 4.15 6.59 8.81 17.09 ▅▇▇▂▁
WS_80mB_mean 0 1 6.69 3.13 0.20 4.21 6.65 8.86 17.20 ▅▇▇▂▁
WS_60m_mean 0 1 6.21 2.94 0.23 3.95 5.96 8.36 16.11 ▃▇▇▂▁
WS_40m_mean 0 1 5.57 2.78 0.20 3.46 5.10 7.68 15.19 ▅▇▆▂▁
WS_20m_mean 0 1 4.68 2.56 0.20 2.59 4.18 6.70 13.27 ▆▇▆▂▁
WD_78m_mean 0 1 230.41 98.56 0.35 136.30 274.10 316.30 359.90 ▁▃▂▂▇
WD_58m_mean 0 1 224.72 99.56 0.03 128.30 269.30 312.30 359.90 ▁▅▂▃▇
temp_80m_mean 0 1 22.01 4.21 11.52 19.23 22.46 24.86 31.18 ▂▃▇▇▂
temp_40m_mean 0 1 21.75 4.47 10.88 18.67 22.07 24.94 31.48 ▂▅▇▆▂
temp_40m_max 0 1 21.89 4.51 11.01 18.79 22.19 25.11 31.76 ▂▅▇▆▂
temp_15m_mean 0 1 21.38 4.72 9.85 17.92 21.63 24.92 31.59 ▂▅▇▇▃
RH_15m_mean 0 1 77.46 17.10 26.87 64.46 81.50 91.80 101.00 ▁▂▃▅▇
P_15m_mean 0 1 1003.77 3.33 996.00 1002.00 1004.00 1006.00 1014.00 ▂▇▆▅▁
DENSIDAD_15m_mean 0 1 30.62 95.98 0.65 1.18 1.23 2.45 350.90 ▇▁▁▁▁

De aquí podemos observar que el conjunto de datos está completo de modo que no hay datos faltantes o NANś que nos puedan causar problemas más adelante.

library(naniar)
## 
## Attaching package: 'naniar'
## The following object is masked from 'package:skimr':
## 
##     n_complete
vis_miss(data)

HABLEMOS DE VELOCIDADES DE VIENTO

Veamos la curva de potencias del aerogenerador que nos interesaría poner en el citio, un modelo Vestas de 2MW. De esta tabla podemos ver que la velocidad de arranque, en la cuál se comienza a producir potencia, es a los 3 m/s. Potencia nominal: 2,000.0 kW Rangos de potencia flexibles: - Velocidad del viento: 4.0 m/s Velocidad nominal del viento: 15.0 m/s Velocidad del viento de corte: 25.0 m/s Velocidad del viento de supervivencia: 60.0 m/s Wind zone (DIBt): - Wind class (IEC): -

cpot<- read.csv("./Datos/Vestas.csv", header=T);cpot
##    Vel..m.s. Pow.kW.
## 1          0       0
## 2          1       0
## 3          2       0
## 4          3      61
## 5          4     156
## 6          5     283
## 7          6     452
## 8          7     669
## 9          8     951
## 10         9    1289
## 11        10    1576
## 12        11    1793
## 13        12    1945
## 14        13    2000
## 15        14    2000
## 16        15    2000
## 17        16    2000
## 18        17    2000
## 19        18    2000
## 20        19    2000
## 21        20    2000
## 22        21    2000
## 23        22    2000
## 24        23    2000
## 25        24    2000
## 26        25       0
## 27        26       0
## 28        27       0
## 29        28       0

Me atrevo a unir los puntos de relaciones de velocidades con potencia generada porque en teoría esta debería ser una curva continua. Ya que para toda velocidad de viento, variable continua, existiría una potencia generada.

plot(cpot[,1], cpot[,2], type = "l",
     col = "red", ylab = "Potencia kW", xlab = "Velocidad de viento m/s")     # Color de la línea

# Agregando puntos
points(cpot[,1], cpot[,2],          # Coordenadas
       pch = 21,      # Símbolo
       cex = 2,       # Tamaño del símbolo
       bg = "green",  # Color de fondo del símbolo
       col = "blue",  # Color del borde del símbolo
       lwd = 3)       # Ancho del borde del símbolo

Podemos decir que la velocidad de corte es a los 25 m/s. Y la potencia nominal es de 2MWatts, se llega a esta potencia a los aproximadamente en el rango [13,25) m/s. Este es el intervalo que nos interesa para que se pueda extraer la mayor energía del viento. Veamos si las velocidades del citio coinciden para que se instale un modelo Vestas 2MW.] se tiene una frecuencia acumulada del 98.06%, de modo que tan solo el 1.94% del tiempo se estuvo en la potencia nominal.

Descartaré las alturas de 20,, ya que este modelo tiene un radio de 40m en sus aspas, lo que interpretamos como la altura mínima a la cual se debe ubicar el aerogenerador.

library(fdth)
(tabla_contm<- fdt(na.omit(WS_60m_mean),breaks="FD",right=T))
##       Class limits   f   rf rf(%)   cf  cf(%)
##  (0.22653,0.76142]  26 0.01  0.60   26   0.60
##   (0.76142,1.2963]  76 0.02  1.76  102   2.36
##    (1.2963,1.8312] 129 0.03  2.99  231   5.35
##    (1.8312,2.3661] 169 0.04  3.91  400   9.26
##     (2.3661,2.901] 186 0.04  4.30  586  13.56
##     (2.901,3.4358] 220 0.05  5.09  806  18.65
##    (3.4358,3.9707] 285 0.07  6.60 1091  25.25
##    (3.9707,4.5056] 294 0.07  6.80 1385  32.05
##    (4.5056,5.0405] 334 0.08  7.73 1719  39.78
##    (5.0405,5.5754] 275 0.06  6.36 1994  46.15
##    (5.5754,6.1103] 240 0.06  5.55 2234  51.70
##    (6.1103,6.6452] 207 0.05  4.79 2441  56.49
##    (6.6452,7.1801] 229 0.05  5.30 2670  61.79
##    (7.1801,7.7149] 239 0.06  5.53 2909  67.32
##    (7.7149,8.2498] 271 0.06  6.27 3180  73.59
##    (8.2498,8.7847] 301 0.07  6.97 3481  80.56
##    (8.7847,9.3196] 237 0.05  5.48 3718  86.04
##    (9.3196,9.8545] 163 0.04  3.77 3881  89.82
##    (9.8545,10.389] 134 0.03  3.10 4015  92.92
##    (10.389,10.924]  85 0.02  1.97 4100  94.89
##    (10.924,11.459]  41 0.01  0.95 4141  95.83
##    (11.459,11.994]  26 0.01  0.60 4167  96.44
##    (11.994,12.529]  30 0.01  0.69 4197  97.13
##    (12.529,13.064]  28 0.01  0.65 4225  97.78
##    (13.064,13.599]  32 0.01  0.74 4257  98.52
##    (13.599,14.134]  23 0.01  0.53 4280  99.05
##    (14.134,14.668]  14 0.00  0.32 4294  99.38
##    (14.668,15.203]  14 0.00  0.32 4308  99.70
##    (15.203,15.738]  10 0.00  0.23 4318  99.93
##    (15.738,16.273]   3 0.00  0.07 4321 100.00

De esto,del sexto renglón, podemos decir que para 60m, el 13.56% de los días tienen velocidades tan vajas que no será suficiente para el arranque, sin embargo, el resto, el 86.44% de los días estará produciendo potencia, pues ni siquiera se detendrá por una velocidad de corte, velocidad en que la maquína se debe detener por seguridad. Sin embargo hasta el intervalo (12.529,13.064] se tiene una frecuencia acumulada del 97.78%, de modo que tan solo el 2.22% del tiempo se estuvo en la potencia nominal.

library(fdth)
(tabla_contm<- fdt(na.omit(WS_40m_mean),breaks="FD",right=T))
##     Class limits   f   rf rf(%)   cf  cf(%)
##  (0.1961,0.7184]  30 0.01  0.69   30   0.69
##  (0.7184,1.2407]  77 0.02  1.78  107   2.48
##   (1.2407,1.763] 158 0.04  3.66  265   6.13
##   (1.763,2.2853] 195 0.05  4.51  460  10.65
##  (2.2853,2.8076] 231 0.05  5.35  691  15.99
##  (2.8076,3.3299] 295 0.07  6.83  986  22.82
##  (3.3299,3.8522] 356 0.08  8.24 1342  31.06
##  (3.8522,4.3745] 361 0.08  8.35 1703  39.41
##  (4.3745,4.8968] 346 0.08  8.01 2049  47.42
##  (4.8968,5.4191] 277 0.06  6.41 2326  53.83
##  (5.4191,5.9414] 218 0.05  5.05 2544  58.88
##  (5.9414,6.4637] 207 0.05  4.79 2751  63.67
##   (6.4637,6.986] 226 0.05  5.23 2977  68.90
##   (6.986,7.5084] 189 0.04  4.37 3166  73.27
##  (7.5084,8.0307] 287 0.07  6.64 3453  79.91
##   (8.0307,8.553] 246 0.06  5.69 3699  85.61
##   (8.553,9.0753] 173 0.04  4.00 3872  89.61
##  (9.0753,9.5976] 126 0.03  2.92 3998  92.52
##   (9.5976,10.12]  92 0.02  2.13 4090  94.65
##   (10.12,10.642]  46 0.01  1.06 4136  95.72
##  (10.642,11.164]  32 0.01  0.74 4168  96.46
##  (11.164,11.687]  28 0.01  0.65 4196  97.11
##  (11.687,12.209]  28 0.01  0.65 4224  97.76
##  (12.209,12.731]  21 0.00  0.49 4245  98.24
##  (12.731,13.254]  25 0.01  0.58 4270  98.82
##  (13.254,13.776]  19 0.00  0.44 4289  99.26
##  (13.776,14.298]  16 0.00  0.37 4305  99.63
##  (14.298,14.821]  12 0.00  0.28 4317  99.91
##  (14.821,15.343]   4 0.00  0.09 4321 100.00

De esto,del sexto renglón, podemos decir que para 40m, el 22.82% de los días tienen velocidades tan vajas que no será suficiente para el arranque, sin embargo, el resto, el 77.18%% de los días estará produciendo potencia, pues ni siquiera se detendrá por una velocidad de corte, velocidad en que la maquína se debe detener por seguridad. Sin embargo hasta el intervalo (12.731,13.254] se tiene una frecuencia acumulada del 98.82%, de modo que tan solo el 1.28% del tiempo se estuvo en la potencia nominal.

En este punto me podría aventurar a decir que el lugar recomendable para poner el aerogenerados es a la altura de 60m, estando en potencia nominal el 2.22% del tiempo. Sin embargo, irémos desmenuzando este punto más adelante.

Hagamos el gráfico de caja y veamos qué podemos decir de la distribución de velocidades de viento.

boxplot(WS_60m_mean,horizontal = T,col="yellow")

quantile(na.omit(WS_60m_mean), probs = c(0.25,0.50,0.75))
##    25%    50%    75% 
## 3.9540 5.9575 8.3592

Para 60m: Podemos decir que la mediana está simétrica al primer y al tercer cuartil. Como el primer cuartil está en 3.95 y el tercero en 8.3, podemos decir que el 50% de las velocidades del viento a 80m está en este rango, suficiente para producir potencia pero insuficiente para la potencia nominal. Las velocidades de viento mayores a 15m/s podemos considerarlas átipicas.

boxplot(WS_80mA_mean,horizontal = T,col="lightgreen")

quantile(na.omit(WS_80mA_mean), probs = c(0.25,0.50,0.75))
##    25%    50%    75% 
## 4.1490 6.5865 8.8064

Para 80m: Podemos decir que la mediana está simétrica al primer y al tercer cuartil. Como el primer cuartil está en 4.14 y el tercero en 8.8, podemos decir que el 50% de las velocidades del viento a 80m está en este rango, suficiente para producir potencia pero insuficiente para la potencia nominal. Las velocidades de viento mayores a 15m/s podemos considerarlas átipicas.

boxplot(WS_40m_mean,horizontal = T,col="red")

quantile(na.omit(WS_40m_mean), probs = c(0.25,0.50,0.75))
##    25%    50%    75% 
## 3.4588 5.1003 7.6850

Del gráfico de cajas podemos cambiar de postura y decir que el 50% de los datos tendrán velocidades más altas a una altura de 80m. Estando en el rango de arranque entre [4.14,8.8] m/s. Pero algo sí que podemos hacer y es descartar la altura de 40m.

Para 80m:

El promedio es 6.58 m/s.

Varianza y desciación estándar:

var(na.omit(WS_80mA_mean)); sd(na.omit(WS_80mA_mean))
## [1] 9.721791
## [1] 3.117979

Coeficiente de asimetría:

(j3<- moments::skewness(WS_80mA_mean,na.rm=T))
## [1] 0.3050912

Coeficiente de curtosis:

(j4<- moments::kurtosis(WS_80mA_mean,na.rm=T))
## [1] 2.733442

Como es menor a tres, podemos decir que será una distribución platycúrtica, o ancha respectro a la distribición normal.

library(BHH2)

dotPlot(na.omit(WS_80mA_mean),main="Diagrama de puntos (muestra aleatoria)",xlab="Velocidad de viento",pch=16)

library(agricolae)
hist_80<- hist(WS_80mA_mean)

ogive.freq(hist_80,main="Ojiva (Velocidades a 80m)",xlab="Límite superior de
           la clase",ylab="Pi (%)")  

##    Límite superior de\n           la clase    RCF
## 1                                        0 0.0000
## 2                                        1 0.0118
## 3                                        2 0.0583
## 4                                        3 0.1282
## 5                                        4 0.2333
## 6                                        5 0.3370
## 7                                        6 0.4464
## 8                                        7 0.5378
## 9                                        8 0.6411
## 10                                       9 0.7707
## 11                                      10 0.8646
## 12                                      11 0.9320
## 13                                      12 0.9586
## 14                                      13 0.9671
## 15                                      14 0.9833
## 16                                      15 0.9912
## 17                                      16 0.9968
## 18                                      17 0.9995
## 19                                      18 1.0000
## 20                                      19 1.0000

Y el diagrama de tallo y ojas. Evaluando hasta el primer decimal:

stem(WS_80mA_mean, scale=1)  # valor por default
## 
##   The decimal point is at the |
## 
##    0 | 2222222223333344455555667777788888888889999999
##    1 | 00000000011111111112222222222222333333333333333333344444444444444444+116
##    2 | 00000000000000000000000000011111111111111111111111111111111111111222+215
##    3 | 00000000000000000000000000000001111111111111111111111111111111111112+356
##    4 | 00000000000000000000000000000000000000000000000000000000000000011111+382
##    5 | 00000000000000000000000000000000000000000111111111111111111111111111+394
##    6 | 00000000000000000000000000000000000111111111111111111111111111111111+314
##    7 | 00000000000000000000000000000000000000001111111111111111111111111111+367
##    8 | 00000000000000000000000000000000000000000111111111111111111111111111+486
##    9 | 00000000000000000000000000000000111111111111111111111111111111111111+317
##   10 | 00000000000000000000000000000000000000001111111111111111111111112222+217
##   11 | 00000000000000000000000000111111111111112222222222222233333333333444+46
##   12 | 0000000011111223334445555555666667778899
##   13 | 00000001111112222233333333444444455566666677777777777778888889999999
##   14 | 000011111111122233334444456667778888
##   15 | 000111123344444457799
##   16 | 000001112245788
##   17 | 01

Se ve claramente que la distribución es asimétrica a la derecha.

Esto lo confirmamos del histograma:

par(mfrow=c(1,2))
hist(WS_80mA_mean,col=rainbow(6), main="Vel. viento a 80m", prob=T)
hist(WS_80mA_mean,col=rainbow(6), main="Vel. viento a 80m", prob=F)

par(mfrow=c(1,2))
hist(WS_60m_mean,col=rainbow(6), main="Vel. viento a 60m", prob=T)
hist(WS_60m_mean,col=rainbow(6), main="Vel. viento a 60m", prob=F)

Variable cuantitativa vs variable cuantitativa: ASOCIACION ENTRE VARIABLES de velocidad de viento

vientos=cbind(WS_80mA_mean,WS_60m_mean,WS_40m_mean,WS_20m_mean)
round(cor(vientos),digits = 2) 
##              WS_80mA_mean WS_60m_mean WS_40m_mean WS_20m_mean
## WS_80mA_mean         1.00        0.99        0.95        0.89
## WS_60m_mean          0.99        1.00        0.98        0.93
## WS_40m_mean          0.95        0.98        1.00        0.97
## WS_20m_mean          0.89        0.93        0.97        1.00
pairs(~WS_80mA_mean+WS_60m_mean+WS_40m_mean+WS_20m_mean, data = mtcars)

De la gráfica de disperción y de las correlaciones podemos concluir que las relaciones van disminuyendo conforme la diferencia en altura. Lo cual tiene sentido, ya que los flujos de aire cambian con la altura.

Interpretemos a las velocidades como una valiable cualitativa.

Como nos interesa la clasificación del viento en sus velocidades de arranque, nominal y de corte, convirtamos este conjunto de datos en una variable cualitativa.

TipoViento80 <- NA
TipoViento80[WS_80mA_mean >= 0 & WS_80mA_mean <= 3] <- "Nulo"   
TipoViento80[WS_80mA_mean> 3 & WS_80mA_mean <= 13] <- "Arranque" 
TipoViento80[WS_80mA_mean> 13 & WS_80mA_mean <= 25] <- "Nominal"
TipoViento80[WS_80mA_mean> 25 & WS_80mA_mean <= 60] <- "Corte"
TipoViento80[WS_80mA_mean> 60] <- "Peligro"



TipoViento60 <- NA
TipoViento60[WS_60m_mean >= 0 & WS_60m_mean <= 3] <- "Nulo"   
TipoViento60[WS_60m_mean> 3 & WS_60m_mean <= 13] <- "Arranque" 
TipoViento60[WS_60m_mean> 13 & WS_60m_mean <= 25] <- "Nominal"
TipoViento60[WS_60m_mean> 25 & WS_60m_mean <= 60] <- "Corte"
TipoViento60[WS_60m_mean> 60] <- "Peligro"
 

TipoViento40 <- NA
TipoViento40[WS_40m_mean >= 0 & WS_40m_mean <= 3] <- "Nulo"   
TipoViento40[WS_40m_mean> 3 & WS_40m_mean <= 13] <- "Arranque" 
TipoViento40[WS_40m_mean> 13 & WS_40m_mean <= 25] <- "Nominal"
TipoViento40[WS_40m_mean> 25 & WS_40m_mean <= 60] <- "Corte"
TipoViento40[WS_40m_mean> 60] <- "Peligro"




TipoViento20 <- NA
TipoViento20[WS_20m_mean >= 0 & WS_20m_mean <= 3] <- "Nulo"   
TipoViento20[WS_20m_mean> 3 & WS_20m_mean <= 13] <- "Arranque" 
TipoViento20[WS_20m_mean> 13 & WS_20m_mean <= 25] <- "Nominal"
TipoViento20[WS_20m_mean> 25 & WS_20m_mean <= 60] <- "Corte"
TipoViento20[WS_20m_mean> 60] <- "Peligro"

Ahora el análicis se centrará en la comparación de variables cualitativas y cuantitativas.

A 80m de altura

# Otra manera
par(mfrow=c(1,1))   # ventana original
library(paletteer)
library(ggthemes)

boxplot(TipoViento80=="Nulo",
        TipoViento80=="Arranque",
        TipoViento80=="Nominal",
        TipoViento80=="Corte",
        TipoViento80=="Peligro",
        col=paletteer_c("ggthemes::Sunset-Sunrise Diverging", 3),names=c("Nulo","Arranque","Nominal","Corte","Peligro"))

De los gráficos de cajas podemos observar que la mayoría de las velocidades de viento se clasifican como velocidades de arranque, es decir, que van de [3,13)m/s. Pasa algo bien curioso y es que en estas clasificaciones la media, el primer cuartil y el tercero parecen coincidir. También podemos decir que las velocidades no son lo suficientemente altas como para que se detenga por seguridad. Y que este parece ser un buen citio si estamos evitando que el aerogenerador se rompa al haber flujos que exedan los 60 m/s.

library(ggplot2)
ggplot(data, 
       aes(x = as.factor(TipoViento80), 
           y = WD_78m_mean)) +
  geom_boxplot() +
 labs(title = "Distribución de la potencia generada por tipo de viento")

Interpretar este gráfici podría ser un reto, debido a que no conocemos a fondo la manera del registro de datos. Mi explicación para los datos nulos que en la gráfica queneran potencia es que al ser tipo arranque la clasificación más común durante un lapso se mantuvo en movimiento las aspas generando potencia, pero solo por un breve lapso de tiempo la velocidad de viento dismunuyó considerablemente a tipo Nulo, sin embargo las aspas ya tenían suficiente momento de inercia, de modo que al momento de registrar la velocidad de viento, las aspas no estaban detenidas del todo por eso se registró potencia. Vamos a continuación cómo esta situación puede ser posible.

head(TipoViento80,50)
##  [1] "Arranque" "Arranque" "Arranque" "Arranque" "Arranque" "Arranque"
##  [7] "Arranque" "Arranque" "Arranque" "Arranque" "Arranque" "Arranque"
## [13] "Arranque" "Arranque" "Arranque" "Arranque" "Arranque" "Arranque"
## [19] "Arranque" "Arranque" "Arranque" "Arranque" "Arranque" "Arranque"
## [25] "Arranque" "Arranque" "Arranque" "Arranque" "Arranque" "Arranque"
## [31] "Arranque" "Arranque" "Arranque" "Arranque" "Arranque" "Arranque"
## [37] "Arranque" "Nulo"     "Arranque" "Arranque" "Arranque" "Nulo"    
## [43] "Arranque" "Arranque" "Arranque" "Arranque" "Arranque" "Arranque"
## [49] "Arranque" "Arranque"
library(ggplot2)

# Distribucion del Ozono 

nivel.viento<- as.factor(TipoViento80)   # para hacerla categorica

ggplot(data, aes(x = DD, 
                     y = WD_78m_mean, 
                     color=nivel.viento)) +
  geom_point() +
  labs(title = "Distribución de la potencia generada por tipo de viento durante noviembre",ylab="Potencia kW")

DD reprencenta los 30 días de noviembre.

library(scales)
dataV<-cbind(TipoViento80,WD_78m_mean)

ggplot(data, 
       aes(y = factor(TipoViento80), 
           x = WD_78m_mean,color = nivel.viento)) +
  geom_jitter() + 
  labs(title = "Potencia registrada por tipo de viento")+
  labs(y= "Tipo de viento", x = "Potencia kW")+
  theme_minimal()

Creo que en este punto es conveniente analizar si la potencia tiene relación con la densidad del aire como dice la teoria y esta última con la temperatura, con la irradiancia y con la humedad. Si pudieramos demostrarlo, pobablemente estas variables nos sirvan para hacer regresión líneal más adelante.

dataDen=cbind(WS_80mA_mean,WD_78m_mean,DENSIDAD_15m_mean)

round(cor(dataDen),digits = 2) 
##                   WS_80mA_mean WD_78m_mean DENSIDAD_15m_mean
## WS_80mA_mean              1.00        0.15              0.00
## WD_78m_mean               0.15        1.00              0.11
## DENSIDAD_15m_mean         0.00        0.11              1.00
pairs(~WS_80mA_mean+WD_78m_mean+DENSIDAD_15m_mean, data = mtcars)

No hay correlación entre la variable densidad y la velocidad del viento, y la correlación entre potencia y densidad es de apenas 0.11.

dataDen=cbind(WS_80mA_mean,WD_78m_mean,DENSIDAD_15m_mean,RH_15m_mean,P_15m_mean,temp_80m_mean)

round(cor(dataDen),digits = 2) 
##                   WS_80mA_mean WD_78m_mean DENSIDAD_15m_mean RH_15m_mean
## WS_80mA_mean              1.00        0.15              0.00       -0.28
## WD_78m_mean               0.15        1.00              0.11       -0.13
## DENSIDAD_15m_mean         0.00        0.11              1.00        0.12
## RH_15m_mean              -0.28       -0.13              0.12        1.00
## P_15m_mean               -0.18       -0.56             -0.09       -0.03
## temp_80m_mean             0.00        0.57              0.19       -0.28
##                   P_15m_mean temp_80m_mean
## WS_80mA_mean           -0.18          0.00
## WD_78m_mean            -0.56          0.57
## DENSIDAD_15m_mean      -0.09          0.19
## RH_15m_mean            -0.03         -0.28
## P_15m_mean              1.00         -0.36
## temp_80m_mean          -0.36          1.00
pairs(~WS_80mA_mean+WD_78m_mean+DENSIDAD_15m_mean+RH_15m_mean+P_15m_mean+temp_80m_mean, data = mtcars)

De los coeficientes de correlación podemos decir que la potencia tiene mayor correlación con la irradiancia y la temperatura. Y por el coeficiente de -0.28 podemos decir que la humedad relativa influye inversamente en la velocidad de viento.

2 Inferencia Estádistica Clásica

Tomaré como variable cuantitativa a la velocidad de viento a 80m. Cuyo histograma es:

hist(WS_80mA_mean,probability=T,col="grey50",
     main="Concentración de ozono en NY",
     xlab="Ozono (ppm)")  # histograma cuyo eje vertical esta en escala
# de "probabilidades".

lines(density(WS_80mA_mean),lwd=2,col="blue")   # curva de densidad de probabilidad

# Para darnos una idea aproximada del perfil
library(fitdistrplus)
## Loading required package: MASS
## Loading required package: survival
plotdist(WS_80mA_mean, histo = TRUE, demp = TRUE)

Realmente no se me viene ninguna distribución a la mente de tan solo ver el perfil de la distribución así que usaré la gráfica de kullen & Frey. Esta propone distribuciones comparando el coeficiente de asímetria y de kurtósis.

descdist(data = WS_80mA_mean, graph = TRUE)

## summary statistics
## ------
## min:  0.18796   max:  17.089 
## median:  6.5865 
## mean:  6.642483 
## estimated sd:  3.117979 
## estimated skewness:  0.3051972 
## estimated kurtosis:  2.734523

De aquí me atrevo a proponer a la distribución betta y weibull. Incluso me arriesgaría con la distribución gamma. Resulta que en el área de energía eólica para las velocidades de viento se le suele asociar la distribución Weibill, veamos si este será el caso.

y<-WS_80mA_mean
library(univariateML)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
## 
##     select
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tibble)
comparacion_aic <- AIC(
  mlbetapr(y),
  mlexp(y),
  mlinvgamma(y),
  mlgamma(y),
  mllnorm(y),
  mlrayleigh(y),
  mlinvgauss(y),
  mlweibull(y),
  mlinvweibull(y)
  #mlbeta(y)
  #mllgamma(y) 

)
comparacion_aic %>% rownames_to_column(var = "distribución") %>% 
  arrange(AIC)
##      distribución df      AIC
## 1    mlweibull(y)  2 21892.51
## 2   mlrayleigh(y)  1 21980.19
## 3      mlgamma(y)  2 22203.54
## 4      mllnorm(y)  2 22964.97
## 5   mlinvgauss(y)  2 23656.75
## 6     mlbetapr(y)  2 23707.42
## 7   mlinvgamma(y)  2 24791.21
## 8        mlexp(y)  1 25007.50
## 9 mlinvweibull(y)  2 25630.41

De aquí se escogería el que tenga el coeficiente de AIC más pequeño. Que serían en este orden, Distribución Weibull, rayleigth y gamma.

Para estimar puntualmete con el método de máxima verosimilitud, usaré esta función:

Los parámetros para la distribución gamma son:

### METODO IIa: estimacion por el "metodo de MLE" (funcion "fitdistr" de R)

library(MASS)
viento_gamma<- fitdistr(na.omit(y), densfun="gamma") # lo hace con MLE
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
viento_gamma
##      shape         rate   
##   3.63573882   0.54733402 
##  (0.07490746) (0.01209270)
viento_weibull<- fitdistr(na.omit(y), densfun="weibull") # lo hace con MLE
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
viento_weibull
##      shape        scale   
##   2.25007875   7.49127814 
##  (0.02720684) (0.05321970)
# Usamos las distribuciones como las nombra univariateML:

par(mfrow=c(1,1))
hist(y,main = "Precipitación en México (mm) durante 2020",
     freq = FALSE)

lines(mlgamma(y), lwd = 2, lty = 2, col = "red")
lines(mlweibull(y), lwd = 2, lty = 3, col="green")
lines(density(y),lwd=2,col="blue")   # curva de densidad de probabilidad
legend(x = 12, y = 0.08, legend = c("gamma", "weibull"),
       col = c("red", "green"), lty = 1:3)
rug(y)   # representacion 1D de los datos

Por el criterio de AIC, la mejor distribución a ajustar concuerda con el que se hace comunmente en la práctica y es la distribución de probabilidad weibull.

2.2 Estimación de intervalos de confianza y pruebas de hpótesiss

Nos interesa inferir si realmente las alturas a los 80m son mayores que a los 60 metros. Pero antes que nada, chequemos si estos datos satisfacen uno de los supuestos básicos para la inferencia cclásica parámetrica, que es el de normalidad.

y=WS_80mA_mean
z=WS_60m_mean
boxplot(y,z,horizontal = T,main="Velocidades de viento",
        names=c("80 metros","60 metros"))

Visualmente la diferencia en el promedio no parece ser avismal.

Verifiquemos si normalidad se cumple:

En las pruebas gráficas entre más se acerque a la línea recta es más probable la normalidad.

En las pruebas estadísticas la hipotesis nula es que la distribución sea normal. Generalmente se hacen varias pruebas. Checa p-value. Recuerda que si es mayor a 0.05 en general no se rechaza hiótesis nula H0, en este caso la aceptaríamos por lo que habría normalidad.

Pruebas para y=Velocidad a 80m de altura

qqnorm(y,pch=16)       # grafica QQ (cuantil-cuantil)
qqline(y,lwd=2,col=2)

library(fBasics)
## Loading required package: timeDate
## 
## Attaching package: 'timeDate'
## The following objects are masked from 'package:agricolae':
## 
##     kurtosis, skewness
## Loading required package: timeSeries
lillieTest(y) # prueba de hipotesis de bondad de ajuste
## 
## Title:
##  Lilliefors (KS) Normality Test
## 
## Test Results:
##   STATISTIC:
##     D: 0.0435
##   P VALUE:
##     < 2.2e-16 
## 
## Description:
##  Mon Jan 31 22:01:52 2022 by user: chicarata

La prueba lilitest() parece indicar que no se cumple la normalidad, ya que con p-value<=0.01 en general se rechaza hipótesis nula.

# Otras pruebas para normalidad
shapiro.test(y)
## 
##  Shapiro-Wilk normality test
## 
## data:  y
## W = 0.98486, p-value < 2.2e-16

De nuevo se rechaza normalidad

Transformamos mediante la técnica de Box-Cox para “corregir” la normalidad.

library(DescTools)
lambda.opt<- BoxCoxLambda(y, method= "loglik", lower=-5, upper=5); lambda.opt
## [1] 0.7
library(nortest)
lillie.test((y^(lambda.opt)-1)/(lambda.opt))     # no hay normalidad
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  (y^(lambda.opt) - 1)/(lambda.opt)
## D = 0.04772, p-value < 2.2e-16

Por lo que transformando mediante BoxCox y con el valor óptimo de lambda sigue sin haber normalidad.

Pruebas para z=Velocidad a 60m de altura

qqnorm(z,pch=16)       # grafica QQ (cuantil-cuantil)
qqline(z,lwd=2,col=2)

lillieTest(z) # prueba de hipotesis de bondad de ajuste
## 
## Title:
##  Lilliefors (KS) Normality Test
## 
## Test Results:
##   STATISTIC:
##     D: 0.0551
##   P VALUE:
##     < 2.2e-16 
## 
## Description:
##  Mon Jan 31 22:01:52 2022 by user: chicarata

La prueba lilitest() parece indicar que no se cumple la normalidad, ya que con p-value<=0.01 en general se rechaza hipótesis nula.

# Otras pruebas para normalidad
shapiro.test(z)
## 
##  Shapiro-Wilk normality test
## 
## data:  z
## W = 0.98118, p-value < 2.2e-16

De nuevo se rechaza normalidad Transformamos mediante la técnica de Box-Cox para “corregir” la normalidad

library(DescTools)
lambda.opt<- BoxCoxLambda(y, method= "loglik", lower=-5, upper=5); lambda.opt
## [1] 0.7
library(nortest)
lillie.test((z^(lambda.opt)-1)/(lambda.opt))   
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  (z^(lambda.opt) - 1)/(lambda.opt)
## D = 0.042945, p-value < 2.2e-16

Por lo que transformando mediante BoxCox y con el valor óptimo de lambda sigue sin haber normalidad.

Como no se llegó a satisfacer el supuesto de normalidad usarémos nuestra mejor aproximación para los intervalos de confianza y pruebas de hipótesis.

Muestras independientes Empecemos abordadndo el problema como muestras independientes, pensando en que pueden llegar a ser independientes si los flujos de viento sontotalmente diferentes a 60m que a 80m. Hidrodinámicamente sabemos que esta suposición sería muy poco probable, pero veamos qué podemos obtener bajo este supuesto.

Como lo que queremos es verificar si efectivamente las velocidades de viento son mayores a 80m, proponemos la hiótesis nula y la alternativa como:

Hipótesis nula: La velocidad promedio a 80m es estadísticamente igual que a los 60m de altura.

Hipótesis alternativa: La velocidad promedio a 80m vs. a los 60m son estadísticamente diferentes.

\[H_0: \mu(80m) = \mu(60m)\] \[Ha: \mu(80m)> \mu(60m)\]

O visto de otra forma: \[H_0: \mu(80m) -\mu(60m)=0\] \[Ha: \mu(80m)- \mu(60m)>0\]

Hagamos la prueba de hipotesis PARAMETRICA para medias con varianza poblacional desconocida. SUPUESTO: normalidad y muestras independientes

Como no conocemos la varianza poblacional usaremos entonces el estádistico t, que se distribuye como una t de student.

prueba1<- t.test(y,z,alternative="greater",mu=0,paired=F,var.equal = T,conf.level = 0.997)
prueba1
## 
##  Two Sample t-test
## 
## data:  y and z
## t = 6.6368, df = 8640, p-value = 1.699e-11
## alternative hypothesis: true difference in means is greater than 0
## 99.7 percent confidence interval:
##  0.2536346       Inf
## sample estimates:
## mean of x mean of y 
##  6.642483  6.209569

El intervalo de confianza es [0.2536346, Inf) y no pasa por el cero, indicando que del 99.7% de todos los intervalos posible no incluyen al cero, lo cual está a favor de la hipótesis alternativa. El estádisctico de prueba es t = 6.6368, y el valor p=p-value = 1.699e-11 de modo que podemos rechazar la hipótesis nula. Por lo que sí hay diferencias en las velocidades de viento.

Cambiemos un poco el desarrollo del problema: Prueba de hipotesis PARAMETRICA para medias con varianza poblacional desconocida. SUPUESTO: normalidad y muestras independientes OJO: los diagramas de caja-brazos no muestran una diferencia en variabilidad. Lo verificamos mediante una prueba de hipotesis:

Aqui la hipótesis nula es: H_0=Es que las varianzas son iguales. El estadístico de prueba es F-test y se distribuye bajo Fisher–Snedecor.

vartest1<- var.test(y,z,ratio=1,alternative="two.sided",conf.level = 0.98)
vartest1
## 
##  F test to compare two variances
## 
## data:  y and z
## F = 1.1222, num df = 4320, denom df = 4320, p-value = 0.0001523
## alternative hypothesis: true ratio of variances is not equal to 1
## 98 percent confidence interval:
##  1.045487 1.204522
## sample estimates:
## ratio of variances 
##           1.122191

Al ser el valor p-value = 0.0001523, las varianzas NO son estadisticamente iguales.

De modo que ahora tenemos una prueba de hipotesis PARAMETRICA para medias con varianza poblacional desconocida. SUPUESTO: normalidad,muestras independientes y varianzas distintas

prueba2<- t.test(y,z,alternative="greater",mu=0,paired=F,var.equal=F, conf.level=0.99)
prueba2
## 
##  Welch Two Sample t-test
## 
## data:  y and z
## t = 6.6368, df = 8611.5, p-value = 1.699e-11
## alternative hypothesis: true difference in means is greater than 0
## 99 percent confidence interval:
##  0.2811403       Inf
## sample estimates:
## mean of x mean of y 
##  6.642483  6.209569

De modo que para muestras independientes y varianzas distintas el intervalo de confianza es [ 0.2811403, Inf] y no pasa por el cero, indicando que del 99%% de todos los intervalos posibles no incluyen al cero, lo cual está a favor de la hipótesis alternativa. El estádisctico de prueba es t = 6.6368, y el valor p, p-value = 1.699e-11 de modo que podemos rechazar la hipótesis nula. Por lo que sí hay diferencias en las velocidades de viento.

Muestras dependientes Ahora supongamos que son muestras dependientes, pensando en que el flujo de aire a 60m modifica al flujo de aire a 80m y viceversa, o bien que estamos hablando del mismo flujo que se comporta diferente según la altura.

Hagamos la prueba de hipotesis PARAMETRICA para medias con varianza poblacional desconocida. SUPUESTO: normalidad y muestras dependientes.

Como no conocemos la varianza poblacional usaremos entonces el estádistico t, que se distribuye como una t de student.

# II. muestras dependientes



prueba3<- t.test(y,z,alternative="two.sided",mu=0,paired=T,var.equal=F,conf.level=0.99)
prueba3
## 
##  Paired t-test
## 
## data:  y and z
## t = 54.414, df = 4320, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 99 percent confidence interval:
##  0.4124116 0.4534156
## sample estimates:
## mean of the differences 
##               0.4329136

El intervalo de confianza es [0.4124116 0.4534156] y no pasa por el cero, indicando que del 99% de todos los intervalos posible no incluyen al cero, lo cual está a favor de la hipótesis alternativa. El estádisctico de prueba es t = 54.414, y el valor p, p-value =2.2e-16 de modo que podemos rechazar la hipótesis nula. Por lo que sí hay diferencias en las velocidades de viento.

Discutamos un poco acerca de cómo deberíamos abordar esta comparación a partir de los intervalos de confianza, una cosa sí es un hecho y es que las velocidades a 80m de altura son mayores que a 60m de altura.

El intervalo de confianza para muestras independientes y varianzas distintas es [ 0.2811403, Inf] y para muestras dependientes y con varianzas distintas es de [0.4124116, 0.4534156] amobos con una certeza del 99%. De esto podemos comprobar que al hacer el tratamiento de dependencia el intervalo de confianza se refuce considerablemente, por lo que debemos quedarnos con el desarrollo dependiente.

Pudimos comprobar que efectivamente las velocidades de viento a 60m son estadísticamente mayores que a 60m, y que las varianzas son estadísticamente distintas. Además de que son muestras dependientes.

Solo por fines didácticos asumamos una probabilidad de rechazar una hipótesis correcta del 8%, construyamos las regiones críticas para rechazar la hipótesis nula. Como no conocemos la varianza poblacional usaremos entonces el estádistico t, que se distribuye como una t de student. Entonces prueba de hipotesis PARAMETRICA para medias con varianza poblacional desconocida. SUPUESTO: normalidad y muestras dependientes con una certeza del 92% es:

La región crítica queda como describe la imágen.

texto_alternativo Entonces el estadístico de prueba:

prueba2<- t.test(y,z,alternative="greater",mu=0,paired=F,var.equal=F, conf.level=0.92)
prueba2
## 
##  Welch Two Sample t-test
## 
## data:  y and z
## t = 6.6368, df = 8611.5, p-value = 1.699e-11
## alternative hypothesis: true difference in means is greater than 0
## 92 percent confidence interval:
##  0.3412545       Inf
## sample estimates:
## mean of x mean of y 
##  6.642483  6.209569

Como el estadístico de prueba es t = 6.6368, usando la región crítica, rechazamos hipótesis nula si t no está entre [-1.878,1.878], por lo que se rechaza hipótesis nula usando el estadístico de prueba.

El intervalo de confianza para muestras dependientes y varianzas distintas con una certeza del 92% es [ 0.3412545, Inf] y para muestras dependientes y con varianzas distintas es de [0.4124116, 0.4534156] a una certeza del 99%. De esto podemos comprobar que al hacer el tratamiento de dependencia el intervalo de confianza al 99% se reduce considerablemente, por lo que debemos quedarnos con el desarrollo con certeza del 99%.

Al ser mayor la velocidad del viento a 80m de altura conviene poner a esa altura el aerogenerador.

library(Hmisc)
## Loading required package: lattice
## Loading required package: Formula
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:DescTools':
## 
##     %nin%, Label, Mean, Quantile
## The following objects are masked from 'package:dplyr':
## 
##     src, summarize
## The following objects are masked from 'package:base':
## 
##     format.pval, units
library(GGally)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

Empecemos con un modelo que contemple las variables:

WS_80m_mean: Velocidad del viento a la altura 80 metros.

WD_78m_mean: La potencia producida en kW/hora. temp_80m_mean: la temperatura en grados centígrados

RH_15m_mean: La humedad relativa

P_15m_mean: Representa la irradiancia. La radiación solar medida en cada una de las estaciones meteorológicas es ofrecida en unidades de potencia y está en vatios por metro cuadrado (W/m²)

DENSIDAD_15m_mean: la densidad del viento [Kg/m^3]. El tipo de variable es cuantitativa y su escala es de razon. Variable continua.

Mostremos la matriz de dispersión de todas las variables involucradas y matriz de correlación para ver si el modelo líneal es justificable.

datos <- data.frame(velocidad=WS_80mA_mean,potencia =WD_78m_mean, densidad = DENSIDAD_15m_mean, humedad = RH_15m_mean, irradiancia =  P_15m_mean, tempratura=temp_80m_mean)
#names(datos) # Nombres de las variables (columnas)
#datos=cbind(WD_78m_mean, DENSIDAD_15m_mean, RH_15m_mean, P_15m_mean, temp_80m_mean)
ggpairs(datos, lower = list(continuous = "smooth"))

Una conbinación que nos podría dar problemas más adelante es la de temperatura con potencia, tienen una correlación de 0.569, y podrían no ser independientes. Aunque recordemos que correlación no indica causalidad.

#WS_80mA_mean ~ WD_78m_mean, DENSIDAD_15m_mean, RH_15m_mean, P_15m_mean, temp_80m_mean
model1<- lm( WD_78m_mean~WS_80mA_mean + DENSIDAD_15m_mean+ RH_15m_mean+ P_15m_mean+ temp_80m_mean, data = data)  # con todos los predictores
model1
## 
## Call:
## lm(formula = WD_78m_mean ~ WS_80mA_mean + DENSIDAD_15m_mean + 
##     RH_15m_mean + P_15m_mean + temp_80m_mean, data = data)
## 
## Coefficients:
##       (Intercept)       WS_80mA_mean  DENSIDAD_15m_mean        RH_15m_mean  
##         1.161e+04          2.554e+00         -6.778e-03          1.368e-02  
##        P_15m_mean      temp_80m_mean  
##        -1.157e+01          1.003e+01
summary(model1)
## 
## Call:
## lm(formula = WD_78m_mean ~ WS_80mA_mean + DENSIDAD_15m_mean + 
##     RH_15m_mean + P_15m_mean + temp_80m_mean, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -300.179  -36.148    6.868   48.145  190.615 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        1.161e+04  3.720e+02  31.203  < 2e-16 ***
## WS_80mA_mean       2.554e+00  3.772e-01   6.770 1.46e-11 ***
## DENSIDAD_15m_mean -6.778e-03  1.179e-02  -0.575    0.565    
## RH_15m_mean        1.368e-02  7.236e-02   0.189    0.850    
## P_15m_mean        -1.157e+01  3.660e-01 -31.620  < 2e-16 ***
## temp_80m_mean      1.003e+01  3.044e-01  32.950  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 71.55 on 4315 degrees of freedom
## Multiple R-squared:  0.4737, Adjusted R-squared:  0.4731 
## F-statistic: 776.7 on 5 and 4315 DF,  p-value: < 2.2e-16

Este modelo nos está dando un valor \(R^2\) de 0.4731, lo cual nos dice que el modelo es capaz de explicar el 47.31% de la variabilidad. + Intercepto: En promedio, la potencia registrada sube 11,610 unidades independientemente de nuestras variables + β1: Por cada unidad adicional de la velocidad, la potencia registrada decrece en promedio por 2.55 unidades. + β2: Por cada cambio en la densidad del viento, la potencia registrada decrece en promedio por -6.778e-03 unidades. + β3: Por cada cambio en la humedad relativa la potencia registrada sube en promedio 1.368e-02 unidades. + β4:Por unidad de cambio en la irradiancia la potencia registrada decrece en promedio -1.157e+01 + β5:Pr unidad de cambio en la temperatura la potencia registrada sube en promedio 10 unidades.

Vemos que lo que menos contribuye a la regresión es la humedad densidad del viento así como la humedad relativa, así que podríamos considerar quitar estas variables para un segundo modelo.

Estudiemos la distribución de probabilidad de cada variable involucrada. NOTA: La línea roja es la curva normal teórica y la línea azul es la densidad empírica.

library(psych)
## 
## Attaching package: 'psych'
## The following object is masked from 'package:Hmisc':
## 
##     describe
## The following objects are masked from 'package:DescTools':
## 
##     AUC, ICC, SD
## The following object is masked from 'package:fBasics':
## 
##     tr
## The following object is masked from 'package:timeSeries':
## 
##     outlier
## The following objects are masked from 'package:scales':
## 
##     alpha, rescale
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
multi.hist(x = datos, dcol = c("blue", "red"), dlty = c("dotted", "solid"),lwd=c(2,1),
           main = c("IQ","Cerebro","Altura","Peso"))

Nos interesa normalidad porque uno de los supuestos es que los errores siguen una distribución normal, y esto se hereda en parte de las varibles regresoras.

Primero haremos un diagrama de caja y brazos para los residuos del modelo.

boxplot(model1[['residuals']],main='Residuos',ylab='Valor residual')

Recordemos que suponemos que los errores de nuestro modelo tienen una distribución normal y en consecuencia sus residuos también, por lo que el diagrama de caja y brazos nos debería mostrar una caja con brazos simétricos, 1er y 3er cuartiles a la misma distancia del cero y una mediana igual a 0.

Verificación de supuestos
plot(model1, which = 1)

En esta gráfica buscamos verificar la normalidad en los residuos. La línea punteada nos indica la línea teórica de una distribución normal y entre mejor nuestros puntos se distribuyan sobre ésta recta, podemos decir que nuestros residuos tienen una distribución normal. Recordemos que ésta gráfica nos perimte revisar si los residuos tienen algún patrón no lineal. Si hay algúna relación no lineal entre las variables regresoras y la respuesta, se verá aquí. En ésta gráfica podemos ver que nuestros residuos muestran un patrón lineal, e inclusive la linea roja se acerca muy bien a la teórica. Por otro lado, notemos que las observaciones 3011, 37 y 3010 se hacen presentes y pueden ser outliers, así que hay que tenerlas en cuenta para las siguientes gráficas.

plot(model1, which = 2)

En esencia tenemos algunos puntos que se distribuyen muy bien sobre la recta, sin embargo en los extremos tenemos valores que se alejan mucho. Esto podría ser un indicio que uno de los supuestos de nuestro modelo no se están cumpliendo, lo cual nos podría indicar que el modelo es incorrecto o la población de donde se extrajeron los datos no es normal. Sin embargo, a consideración personal, olvidándonos de las observaciones más alejadas, nuestros residuos parecen tener una distribución adecuada. De nuevo las observaciones 37 y 3011,3010 se hacen presentes.

plot(model1, which = 3)

Con ésta gráfica buscamos verificar la homocedasticidad. Ésta gráfica nos permite verificar el supuesto de varianza constante en los resiudos. Como podemos ver,la línea roja tiene un comportamiento casi teórico, pues esta linea idealmente debe de ser horizontal. Además por otro lado, nuestros puntos se presetan como una nube en desorden pero no cubren todo el espacio de forma uniforme. Por lo que podemos concluir que nuestro modelo efectivamente tiene varianza constante en los residuos. Igualmente las observaciones 3011, 37 y 3010 están presentes como observaciones que pueden estar afectando el modelo.

plot(model1, which = 5)

Ésta gráfica nos sirve para determinar si hay datos influyentes en la regresión. Graficamente no puedo dar una conclución.

Elijamos un segundo modelo, según el criterio de Akaike tenemos que las mejores variables para nuestro modelo son:

step(object = model1, direction = "both", trace = 1)
## Start:  AIC=36910.5
## WD_78m_mean ~ WS_80mA_mean + DENSIDAD_15m_mean + RH_15m_mean + 
##     P_15m_mean + temp_80m_mean
## 
##                     Df Sum of Sq      RSS   AIC
## - RH_15m_mean        1       183 22089044 36909
## - DENSIDAD_15m_mean  1      1693 22090554 36909
## <none>                           22088862 36910
## - WS_80mA_mean       1    234621 22323483 36954
## - P_15m_mean         1   5118258 27207120 37809
## - temp_80m_mean      1   5557931 27646793 37878
## 
## Step:  AIC=36908.53
## WD_78m_mean ~ WS_80mA_mean + DENSIDAD_15m_mean + P_15m_mean + 
##     temp_80m_mean
## 
##                     Df Sum of Sq      RSS   AIC
## - DENSIDAD_15m_mean  1      1542 22090586 36907
## <none>                           22089044 36909
## + RH_15m_mean        1       183 22088862 36910
## - WS_80mA_mean       1    259063 22348107 36957
## - P_15m_mean         1   5375169 27464213 37848
## - temp_80m_mean      1   6450363 28539407 38014
## 
## Step:  AIC=36906.83
## WD_78m_mean ~ WS_80mA_mean + P_15m_mean + temp_80m_mean
## 
##                     Df Sum of Sq      RSS   AIC
## <none>                           22090586 36907
## + DENSIDAD_15m_mean  1      1542 22089044 36909
## + RH_15m_mean        1        32 22090554 36909
## - WS_80mA_mean       1    259053 22349639 36955
## - P_15m_mean         1   5373631 27464217 37846
## - temp_80m_mean      1   6615076 28705662 38037
## 
## Call:
## lm(formula = WD_78m_mean ~ WS_80mA_mean + P_15m_mean + temp_80m_mean, 
##     data = data)
## 
## Coefficients:
##   (Intercept)   WS_80mA_mean     P_15m_mean  temp_80m_mean  
##     11621.727          2.530        -11.584          9.981

WS_80m_mean: Velocidad del viento a la altura 80 metros.

WD_78m_mean: La potencia producida en kW/hora. temp_80m_mean: la temperatura en grados centígrados

P_15m_mean: Representa la irradiancia. La radiación solar medida en cada una de las estaciones meteorológicas es ofrecida en unidades de potencia y está en vatios por metro cuadrado (W/m²)

Como habíamos predicho anteriormente, se sugiere quitar a la humedad relativa y la densidad de viento como variables regresoras.

Entonces analicemos este modelo

#WS_80mA_mean ~ WD_78m_mean, DENSIDAD_15m_mean, RH_15m_mean, P_15m_mean, temp_80m_mean
model2<- lm( WD_78m_mean~WS_80mA_mean +  P_15m_mean+ temp_80m_mean, data = data)  # con todos los predictores
summary(model2)
## 
## Call:
## lm(formula = WD_78m_mean ~ WS_80mA_mean + P_15m_mean + temp_80m_mean, 
##     data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -299.886  -36.252    6.919   48.149  190.312 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   11621.7266   361.5677  32.143  < 2e-16 ***
## WS_80mA_mean      2.5301     0.3556   7.115  1.3e-12 ***
## P_15m_mean      -11.5842     0.3575 -32.406  < 2e-16 ***
## temp_80m_mean     9.9812     0.2776  35.955  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 71.53 on 4317 degrees of freedom
## Multiple R-squared:  0.4736, Adjusted R-squared:  0.4733 
## F-statistic:  1295 on 3 and 4317 DF,  p-value: < 2.2e-16

Cemos que con este modelo se pueden describir mayor parte de la muestra, un 47.33%.

plot(model2)

Los supuestos parecen ser los mismos que en el modelo 1.

Veamos cuál variable tiene datos atípicos:

X2 <- data.frame(cbind(WD_78m_mean, P_15m_mean, temp_80m_mean));

X2.labels <- c(WD_78m_mean = "Potencia",
           P_15m_mean = "Irradiancia",
           temp_80m_mean = "Temperatura")
label(X2) <- as.list(X2.labels[match(names(X2), names(X2.labels))])



label <- 1
for (variable in X2){
  boxplot(variable, id.method = 'y', ylab = names(X2)[label], main = X2.labels[label])
  label <- label + 1
}

n<- nrow(X2)
dif.betas <- abs(dfbeta(model2)) > sqrt(2/n)

for (i in 1:ncol(dif.betas)){
  if (i != 1){
    influyentes <- which(dif.betas[, i] == TRUE)
  }
  else{
    influyentes <- which(dif.betas[, 1] == FALSE)
  }
  print(X2.labels[i-1])
  print(names(influyentes))
}
## named character(0)
##   [1] "135"  "232"  "233"  "267"  "268"  "269"  "270"  "271"  "272"  "370" 
##  [11] "409"  "458"  "526"  "527"  "571"  "576"  "581"  "584"  "586"  "699" 
##  [21] "823"  "825"  "966"  "980"  "988"  "996"  "999"  "1002" "1103" "1106"
##  [31] "1153" "1154" "1155" "1156" "1157" "1170" "1194" "1195" "1304" "1320"
##  [41] "1350" "1637" "1642" "1689" "1991" "1997" "2091" "2093" "2094" "2137"
##  [51] "2138" "2139" "2140" "2141" "2142" "2144" "2145" "2146" "2147" "2248"
##  [61] "2270" "2788" "2789" "2790" "2791" "3212" "3213" "3270" "3275" "3280"
##  [71] "3281" "3294" "3386" "3649" "3765" "3767" "3768" "3769" "3770" "3771"
##  [81] "3772" "3774" "3775" "3776" "3777" "3778" "3779" "3783" "3784" "3785"
##  [91] "3786" "3792" "3834" "3875" "3973" "3976" "3977" "3978" "3983" "3991"
## [101] "3995" "4134" "4209"
## WD_78m_mean 
##  "Potencia" 
##  [1] "37"   "40"   "41"   "42"   "43"   "44"   "45"   "46"   "47"   "48"  
## [11] "49"   "50"   "51"   "52"   "53"   "54"   "56"   "57"   "58"   "59"  
## [21] "60"   "61"   "62"   "63"   "1057" "1531" "1596" "1609" "1610" "1618"
## [31] "1619" "1620" "1627" "2178" "2215" "2589" "2590" "2591" "2592" "2593"
## [41] "2594" "2595" "3007" "3008" "3010" "3011" "3013" "3014" "3016" "3017"
## [51] "3018" "3019" "3028" "3029" "3030" "3031" "3034" "3035" "3036" "3037"
## [61] "3038" "3039" "3040" "3041" "3042" "3043" "3044" "3045" "3046" "3047"
## [71] "3048" "3049" "3050" "3051" "3052" "3053" "3054" "3055" "3056" "3057"
## [81] "4070" "4181" "4187"
##    P_15m_mean 
## "Irradiancia" 
##  [1] "37"   "40"   "41"   "42"   "43"   "44"   "45"   "46"   "47"   "48"  
## [11] "49"   "50"   "51"   "52"   "53"   "54"   "55"   "56"   "57"   "1803"
## [21] "1804" "3007" "3008" "3010" "3011" "3013" "3014" "3015" "3016" "3017"
## [31] "3018" "3019" "3020" "3021" "3022" "3023" "3344" "3360"
## temp_80m_mean 
## "Temperatura" 
## character(0)
remove <- c( "37" ,  "40" ,  "41" ,  "42"  , "43" ,  "44"  , "45"  , "46" ,  "47"  , "48",   "49","50" ,  "51" ,  "52"  , "53"  , "54"   ,"55"  , "56"  , "57" ,  "1803", "1804", "3007","3008" ,"3010" ,"3011", "3013", "3014" ,"3015", "3016" ,"3017" ,"3018", "3019", "3020","3021" ,"3022", "3023", "3344", "3360","58" ,  "59"  , "3024")

data.nuevo <- data[rownames(data) %nin% remove, ]
attach(data.nuevo)
## The following objects are masked from data:
## 
##     DD, DENSIDAD_15m_mean, hh, mm, MM, P_15m_mean, RH_15m_mean,
##     temp_15m_mean, temp_40m_max, temp_40m_mean, temp_80m_mean,
##     WD_58m_mean, WD_78m_mean, WS_20m_mean, WS_40m_mean, WS_60m_mean,
##     WS_80mA_mean, WS_80mB_mean
X3 <- data.frame(cbind(WD_78m_mean, P_15m_mean, temp_80m_mean));

X3.labels <- c(WD_78m_mean = "Potencia",
           P_15m_mean = "Irradiancia",
           temp_80m_mean = "Temperatura")
label(X3) <- as.list(X3.labels[match(names(X3), names(X3.labels))])
model3<- lm( WD_78m_mean~WS_80mA_mean +  P_15m_mean+ temp_80m_mean, data = data.nuevo)  # con todos los predictores
summary(model3)
## 
## Call:
## lm(formula = WD_78m_mean ~ WS_80mA_mean + P_15m_mean + temp_80m_mean, 
##     data = data.nuevo)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -268.16  -35.04    5.73   46.34  184.86 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   12992.6743   346.0859  37.542  < 2e-16 ***
## WS_80mA_mean      1.5431     0.3375   4.573 4.95e-06 ***
## P_15m_mean      -12.9331     0.3421 -37.802  < 2e-16 ***
## temp_80m_mean     9.6180     0.2622  36.689  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 67.36 on 4276 degrees of freedom
## Multiple R-squared:  0.5176, Adjusted R-squared:  0.5172 
## F-statistic:  1529 on 3 and 4276 DF,  p-value: < 2.2e-16

!Vaya! Ahora podemos describir a un 51.72% de los datos.

plot(model3)

n<- nrow(X3)
dif.betas <- abs(dfbeta(model3)) > sqrt(2/n)

for (i in 1:ncol(dif.betas)){
  if (i != 1){
    influyentes <- which(dif.betas[, i] == TRUE)
  }
  else{
    influyentes <- which(dif.betas[, 1] == FALSE)
  }
  print(X3.labels[i-1])
  print(names(influyentes))
}
## named character(0)
##  [1] "23"   "139"  "233"  "238"  "267"  "268"  "269"  "270"  "271"  "272" 
## [11] "374"  "529"  "583"  "656"  "823"  "825"  "865"  "946"  "966"  "989" 
## [21] "994"  "995"  "997"  "998"  "999"  "1000" "1001" "1101" "1127" "1153"
## [31] "1154" "1155" "1156" "1157" "1170" "1194" "1195" "1320" "1981" "1992"
## [41] "1999" "2001" "2091" "2093" "2094" "2096" "2137" "2138" "2139" "2140"
## [51] "2141" "2142" "2143" "2144" "2145" "2146" "2147" "2270" "2271" "2469"
## [61] "2624" "2789" "2790" "2791" "3212" "3213" "3269" "3271" "3278" "3293"
## [71] "3402" "3772" "3773" "3776" "3791" "3796" "3874" "3970" "3971" "3977"
## [81] "3980" "3983" "3991" "3995" "3996" "4274"
## WD_78m_mean 
##  "Potencia" 
##  [1] "60"   "61"   "62"   "63"   "1057" "1531" "1582" "1596" "1597" "1601"
## [11] "1602" "1603" "1605" "1608" "1609" "1610" "1618" "1619" "1620" "1622"
## [21] "1627" "2178" "2215" "2589" "2590" "2591" "2592" "2593" "2594" "3028"
## [31] "3029" "3030" "3031" "3034" "3035" "3036" "3037" "3038" "3039" "3040"
## [41] "3041" "3042" "3043" "3044" "3045" "3046" "3047" "3048" "3049" "3050"
## [51] "3051" "3052" "3053" "3055" "3056" "3057" "4070" "4181" "4187" "4194"
## [61] "4195"
##    P_15m_mean 
## "Irradiancia" 
## [1] "2904"
## temp_80m_mean 
## "Temperatura" 
## character(0)
remove <- c( "60" ,  "61"  , "62" ,  "63" ,  "1057", "1531" ,"1582" ,"1596" ,"1597" ,"1601" ,"1602","1603", "1605" ,"1608" ,"1609", "1610", "1618" ,"1619", "1620", "1622", "1627" ,"2178","2215", "2589" ,"2590", "2591", "2592" ,"2593", "2594", "3028", "3029", "3030" ,"3031", "3034", "3035" ,"3036" ,"3037" ,"3038", "3039" ,"3040","3041", "3042", "3043" ,"3044", "3045" ,"3046" ,"3047" ,"3048", "3049" ,"3050" ,"3051" ,"3052" ,"3053", "3055", "3056", "3057" ,"4070" ,"4181", "4187" ,"4194" ,"4195")

data.nuevo <- data[rownames(data) %nin% remove, ]
attach(data.nuevo)
## The following objects are masked from data.nuevo (pos = 3):
## 
##     DD, DENSIDAD_15m_mean, hh, mm, MM, P_15m_mean, RH_15m_mean,
##     temp_15m_mean, temp_40m_max, temp_40m_mean, temp_80m_mean,
##     WD_58m_mean, WD_78m_mean, WS_20m_mean, WS_40m_mean, WS_60m_mean,
##     WS_80mA_mean, WS_80mB_mean
## The following objects are masked from data:
## 
##     DD, DENSIDAD_15m_mean, hh, mm, MM, P_15m_mean, RH_15m_mean,
##     temp_15m_mean, temp_40m_max, temp_40m_mean, temp_80m_mean,
##     WD_58m_mean, WD_78m_mean, WS_20m_mean, WS_40m_mean, WS_60m_mean,
##     WS_80mA_mean, WS_80mB_mean
X4 <- data.frame(cbind(WD_78m_mean, P_15m_mean, temp_80m_mean));

X4.labels <- c(WD_78m_mean = "Potencia",
           P_15m_mean = "Irradiancia",
           temp_80m_mean = "Temperatura")
label(X4) <- as.list(X4.labels[match(names(X4), names(X4.labels))])
model4<- lm( WD_78m_mean~WS_80mA_mean +  P_15m_mean+ temp_80m_mean, data = data.nuevo)  # con todos los predictores
summary(model4)
## 
## Call:
## lm(formula = WD_78m_mean ~ WS_80mA_mean + P_15m_mean + temp_80m_mean, 
##     data = data.nuevo)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -302.990  -36.135    5.793   46.767  187.630 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   11931.7177   350.6936  34.023  < 2e-16 ***
## WS_80mA_mean      2.8961     0.3574   8.104 6.87e-16 ***
## P_15m_mean      -11.8893     0.3467 -34.292  < 2e-16 ***
## temp_80m_mean     9.8042     0.2690  36.450  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 69.17 on 4256 degrees of freedom
## Multiple R-squared:  0.4964, Adjusted R-squared:  0.4961 
## F-statistic:  1399 on 3 and 4256 DF,  p-value: < 2.2e-16

Acá notamos que bajó el porcentaje de datos descritos a un 49.61%

plot(model4)

Los supuestos parecen mantenerce en este modelo.

n<- nrow(X2)
dif.betas <- abs(dfbeta(model2)) > sqrt(2/n)

for (i in 1:ncol(dif.betas)){
  if (i != 1){
    influyentes <- which(dif.betas[, i] == TRUE)
  }
  else{
    influyentes <- which(dif.betas[, 1] == FALSE)
  }
  print(X2.labels[i-1])
  print(names(influyentes))
}
## named character(0)
##   [1] "135"  "232"  "233"  "267"  "268"  "269"  "270"  "271"  "272"  "370" 
##  [11] "409"  "458"  "526"  "527"  "571"  "576"  "581"  "584"  "586"  "699" 
##  [21] "823"  "825"  "966"  "980"  "988"  "996"  "999"  "1002" "1103" "1106"
##  [31] "1153" "1154" "1155" "1156" "1157" "1170" "1194" "1195" "1304" "1320"
##  [41] "1350" "1637" "1642" "1689" "1991" "1997" "2091" "2093" "2094" "2137"
##  [51] "2138" "2139" "2140" "2141" "2142" "2144" "2145" "2146" "2147" "2248"
##  [61] "2270" "2788" "2789" "2790" "2791" "3212" "3213" "3270" "3275" "3280"
##  [71] "3281" "3294" "3386" "3649" "3765" "3767" "3768" "3769" "3770" "3771"
##  [81] "3772" "3774" "3775" "3776" "3777" "3778" "3779" "3783" "3784" "3785"
##  [91] "3786" "3792" "3834" "3875" "3973" "3976" "3977" "3978" "3983" "3991"
## [101] "3995" "4134" "4209"
## WD_78m_mean 
##  "Potencia" 
##  [1] "37"   "40"   "41"   "42"   "43"   "44"   "45"   "46"   "47"   "48"  
## [11] "49"   "50"   "51"   "52"   "53"   "54"   "56"   "57"   "58"   "59"  
## [21] "60"   "61"   "62"   "63"   "1057" "1531" "1596" "1609" "1610" "1618"
## [31] "1619" "1620" "1627" "2178" "2215" "2589" "2590" "2591" "2592" "2593"
## [41] "2594" "2595" "3007" "3008" "3010" "3011" "3013" "3014" "3016" "3017"
## [51] "3018" "3019" "3028" "3029" "3030" "3031" "3034" "3035" "3036" "3037"
## [61] "3038" "3039" "3040" "3041" "3042" "3043" "3044" "3045" "3046" "3047"
## [71] "3048" "3049" "3050" "3051" "3052" "3053" "3054" "3055" "3056" "3057"
## [81] "4070" "4181" "4187"
##    P_15m_mean 
## "Irradiancia" 
##  [1] "37"   "40"   "41"   "42"   "43"   "44"   "45"   "46"   "47"   "48"  
## [11] "49"   "50"   "51"   "52"   "53"   "54"   "55"   "56"   "57"   "1803"
## [21] "1804" "3007" "3008" "3010" "3011" "3013" "3014" "3015" "3016" "3017"
## [31] "3018" "3019" "3020" "3021" "3022" "3023" "3344" "3360"
## temp_80m_mean 
## "Temperatura" 
## character(0)

En un siguiente modelo eliminemos los datos atípicos de la irradiancia que parecen ser los más influyendes del diagrama de caja y brazos.

Ecuación Inferencia individual y conjunta \(R^2\) Supuestos
Potencia+densidad+humedad+irradiancia+temperatura Nintendo 0.4731 Normalidad
Homogeneidad
Potenciairradiancia+temperatura Atlus 0.4733 Normalidad
Homogeneidad
Potenciai+rradiancia+temperatura borrando outliers de irradiancia Atlus 0.5271 Normalidad
Homogeneidad
Potencia+irradiancia+temperatura borrando outliers de potencia Atlus 0.4961 Normalidad
Homogeneidad